
%% plot data and secular variation

S = dir('/Users/manojnair/data/obs_mag_data/BGS data/*.nc');


% selected_obs = ['AAA';	'ABG';	'ABK';	'ASP';	'BDV';	'BEL';	'BFO';	...
%     'BLC';	'BMT';	'BOX';	'BOU'; 'CBB';	'CLF';	'CTA';	'CZT';	'DRV';	...
%     'EBR';	'EYR';	'FCC';	'FUR';	'GNA'; 'GLN';	'HLP';	'HRN';	'KAK';	...
%     'KDU';	'KNY';	'LRM';	'LVV';	'MAW';	'MBO';	'MCQ';	'MEA';	...
%     'MMB';	'NGK';	'NUR';	'NVS';	'OTT';	'PAF';	'RES';	'SOD';	...
%     'STP';	'STJ';	'THY'; 'TUC';	'UPS';	'VIC'];

selected_obs = cell2mat(textread('/Users/manojnair/projects/wmm2015_validation/selected_obs.txt','%s'));

time_array_hour = (datenum(2010,1,1,0,0,30): (1/(24)): datenum(2014,12,31,23,59,30))';

for i = 1:length(S),
    filestring(i,:) = S(i).name(1:3);
end;

[data_obs, ia, ib ] = unique(filestring,'rows');

[spline_obs, ia, ib]= intersect(data_obs, selected_obs, 'rows');

% spline_obs = filestring; % using all observatories

% plot

f = figure(1);
set(f,'Position',[40          19        1557         937]);


annual_time_1 = datenum(2010,1:48,15);
annual_time_2 = datenum(2011,1:48,15);
actual_sec_time = datenum(2010,7:54,17);
annual_time = datenum(2010,1:60,15);

% wmm_fixed_file = '/Users/manojnair/projects/wmm2015_validation/WMM2015_fixed.txt';
% wmm_30days_file = '/Users/manojnair/projects/wmm2015_validation/WMM2015_30days.txt';
%wmm_bgs = '/Users/manojnair/projects/wmm2015_validation/swarm_oer_parent_fixed_200.txt';

wmm_bgs = '/Users/manojnair/projects/wmm2015_validation/WMM2015_coeffs_Version_BGS_20141106.txt';
wmm_ngdc = '/Users/manojnair/projects/wmm2015_validation/WMM2015_NGDC_20141104.txt';

for i = 1:length(spline_obs),
    
    ncfname = ['/Users/manojnair/data/obs_mag_data/BGS data/' spline_obs(i,:) '_1995_2014.nc'];
    
    
    %display(ncfname);
    
    [x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj] = read_geomag_netcdf(ncfname, 131496, 43824, 0);
    
    if any(x_data),
        
        
        
        % fit splines
        L = isnan(x_data);
        
        data_array = x_data(~L);
        first_index = find(~L, 1, 'first');
        
        last_index = find(~L, 1, 'last');
        
        time_array = time_array_hour(~L);
        
       % b1 = min(time_array):365:max(time_array)+10;
       b1 = min(time_array)+(365/2):365:max(time_array)+10;
        
        if length(b1) > 1,
            
            sp_x=spline(b1,x_data(~L)'/spline(b1,eye(length(b1)),time_array'));
            sp_y=spline(b1,y_data(~L)'/spline(b1,eye(length(b1)),time_array'));
            sp_z=spline(b1,z_data(~L)'/spline(b1,eye(length(b1)),time_array'));
            
            if  obj.geospatial_lon > 180,
                obj.geospatial_lon = obj.geospatial_lon - 360;
            end;
            
            
            icount = 1;
            for decyear = 2010:0.1:2014.9,
                [wmm_2010(icount,:)] = wrldmagm(0,obj.geospatial_lat,obj.geospatial_lon,decyear);
                
                geoc_lat = geod2geoc(obj.geospatial_lat, 0, 'WGS84'); % geocentric lat
                geoc_r = ( geocradius(geoc_lat,'WGS84') - 6371.2*1000 )/1000; % geocentric alt in km
                [wmm_2015_ngdc(icount,:)] = magsynth(geoc_r,geoc_lat,obj.geospatial_lon,decyear,wmm_ngdc );
                [wmm_2015_bgs(icount,:)] = magsynth(geoc_r,geoc_lat,obj.geospatial_lon,decyear,wmm_bgs );
                
                icount = icount + 1;
            end;
            
            % secular variation is linear and same through out the model so
            % no need to calculate it through the loop.
            
            
            geoc_lat = geod2geoc(obj.geospatial_lat, 0, 'WGS84'); % geocentric lat
            geoc_r = ( geocradius(geoc_lat,'WGS84') - 6371.2*1000 )/1000; % geocentric alt in km
            
            
            sec_wmm_ngdc = magsynth(geoc_r,geoc_lat, ...
                obj.geospatial_lon,2014.8,wmm_ngdc)-magsynth(geoc_r,geoc_lat, ...
                obj.geospatial_lon,2013.8,wmm_ngdc);
            
            sec_wmm_bgs = magsynth(geoc_r,geoc_lat, ...
                obj.geospatial_lon,2014.8,wmm_bgs)-magsynth(geoc_r,geoc_lat, ...
                obj.geospatial_lon,2013.8,wmm_bgs);
            
            sec_wmm2010 = (wrldmagm(0,obj.geospatial_lat, ...
                obj.geospatial_lon,2014.0)-wrldmagm(0, ...
                obj.geospatial_lat,obj.geospatial_lon,2010))/4;
            
            sec_x_obs =  ppval(sp_x,annual_time_2) - ppval(sp_x,annual_time_1) ;
            sec_y_obs =  ppval(sp_y,annual_time_2) - ppval(sp_y,annual_time_1) ;
            sec_z_obs =  ppval(sp_z,annual_time_2) - ppval(sp_z,annual_time_1);
            
            
            sec_x_obs_eval =  ppval(sp_x,datenum(2014,6,19)) - ppval(sp_x,datenum(2013,6,19)) ;
            sec_y_obs_eval =  ppval(sp_y,datenum(2014,6,19)) - ppval(sp_y,datenum(2013,6,19)) ;
            sec_z_obs_eval =  ppval(sp_z,datenum(2014,6,19)) - ppval(sp_z,datenum(2013,6,19));
            
            fprintf('%s %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f \n', spline_obs(i,:), ...
                obj.geospatial_lat,obj.geospatial_lon,  sec_x_obs_eval - sec_wmm_ngdc(1), ...
                sec_y_obs_eval - sec_wmm_ngdc(2),sec_z_obs_eval - sec_wmm_ngdc(3),...
                sec_x_obs_eval - sec_wmm_bgs(1), ...
                sec_y_obs_eval - sec_wmm_bgs(2),sec_z_obs_eval - sec_wmm_bgs(3));
            
            
            delta_sec(i,1) = obj.geospatial_lat;
            delta_sec(i,2) = obj.geospatial_lon;
            delta_sec(i,3) = sec_x_obs(end) - sec_wmm_bgs(1);
            delta_sec(i,4) = sec_y_obs(end) - sec_wmm_bgs(2);
            delta_sec(i,5) = sec_z_obs(end) - sec_wmm_bgs(3);
            delta_sec(i,6) = sec_x_obs(end) - sec_wmm_ngdc(1);
            delta_sec(i,7) = sec_y_obs(end) - sec_wmm_ngdc(2);
            delta_sec(i,8) = sec_z_obs(end) - sec_wmm_ngdc(3);
            
            
            % remove non data region.
            
%             sec_x_obs(actual_sec_time > time_array_hour(last_index)) = NaN;
%             sec_y_obs(actual_sec_time > time_array_hour(last_index)) = NaN;
%             sec_z_obs(actual_sec_time > time_array_hour(last_index)) = NaN;
%             
%             
%             
%             
%             f = figure(1);
%             set(f,'Position',[ 560    39   778   909]);
%             subplot(321);
%             
%             plot(actual_sec_time, sec_x_obs,'LineWidth',2);
%             set(gca,'FontSize',16);
%             hold on;
%             plot(datenum(2010,1,1), sec_wmm2010(1),'r*','MarkerSIze',20);
%             plot(datenum(2015,1,1), sec_wmm_ngdc(1),'k*','MarkerSIze',20);
%             plot(datenum(2015,1,1), sec_wmm_bgs(1),'b*','MarkerSIze',20);
%             % plot(datenum(2014,7,2), ppval(sp_x,datenum(2014,7,2)),'g*','MarkerSIze',20);
%             
%             
%             h=legend('OBS','WMM2010','WMM2015NGDC','WMM2015BGS');
%             set(h,'FontSize',8,'Location','NorthWest');
%             
%             
%             datetick('x','keeplimits');
%             ylabel('X(nT/yr)');
%             title(sprintf('%s lat %5.2f lon %5.2f ',spline_obs(i,:), ...
%                 obj.geospatial_lat, obj.geospatial_lon));
%             
%             subplot(323);
%             plot(actual_sec_time, sec_y_obs,'LineWidth',2);
%             set(gca,'FontSize',16);
%             hold on;
%             plot(datenum(2010,1,1), sec_wmm2010(2),'r*','MarkerSIze',20);
%             plot(datenum(2015,1,1), sec_wmm_ngdc(2),'k*','MarkerSIze',20);
%             plot(datenum(2015,1,1), sec_wmm_bgs(2),'b*','MarkerSIze',20);
%             %plot(datenum(2014,7,2), ppval(sp_y,datenum(2014,7,2)),'g*','MarkerSIze',20);
%             
%             
%             datetick('x','keeplimits');
%             ylabel('Y(nT/yr)');
%             
%             subplot(325);
%             plot(actual_sec_time, sec_z_obs,'LineWidth',2);
%             set(gca,'FontSize',16);
%             hold on;
%             plot(datenum(2010,1,1), sec_wmm2010(3),'r*','MarkerSIze',20);
%             plot(datenum(2015,1,1), sec_wmm_ngdc(3),'k*','MarkerSIze',20);
%             plot(datenum(2015,1,1), sec_wmm_bgs(3),'b*','MarkerSIze',20);
%             % plot(datenum(2014,7,2), ppval(sp_z,datenum(2014,7,2)),'g*','MarkerSIze',20);
%             
%             
%             datetick('x','keeplimits');
%             ylabel('Z(nT/yr)');
%             
%             
%             
%             subplot(322);
%             
%             plot(time_array_hour,x_data);
%             hold on;
%             spline_fit = ppval(time_array_hour,sp_x);
%             spline_fit(1:first_index) = NaN;
%             spline_fit(last_index:end) = NaN;
%             plot(time_array_hour,spline_fit,'r');
%             plot(datenum(2010:0.1:2014.9,0,0), wmm_2015_ngdc(:,1),'r.-','MarkerSize',15);
%             plot(datenum(2010:0.1:2014.9,0,0), wmm_2015_bgs(:,1),'b.-','MarkerSize',15);
%             
%             h=legend('OBS','FIT','WMM2015NGDC','WMM2015BGS');
%             set(h,'FontSize',8);
%             set(gca,'FontSize',16);
%             
%             ylabel('X nT');
%             axis([datenum(2010,1,1), datenum(2015,2,1), -inf inf]);
%             datetick('x','keeplimits');
%             
%             
%             subplot(324);
%             
%             plot(time_array_hour,y_data);
%             hold on;
%             spline_fit = ppval(time_array_hour,sp_y);
%             spline_fit(1:first_index) = NaN;
%             spline_fit(last_index:end) = NaN;
%             plot(time_array_hour,spline_fit,'r');
%             
%             hold on;
%             plot(datenum(2010:0.1:2014.9,0,0), wmm_2015_ngdc(:,2),'r.-','MarkerSize',15);
%             plot(datenum(2010:0.1:2014.9,0,0), wmm_2015_bgs(:,2),'b.-','MarkerSize',15);
%             
%             
%             set(gca,'FontSize',16);
%             ylabel('Y nT');
%             axis([datenum(2010,1,1), datenum(2015,1,1), -inf inf]);
%             datetick('x','keeplimits');
%             subplot(326);
%             
%             plot(time_array_hour,z_data);
%             hold on;
%             spline_fit = ppval(time_array_hour,sp_z);
%             spline_fit(1:first_index) = NaN;
%             spline_fit(last_index:end) = NaN;
%             plot(time_array_hour,spline_fit,'r');
%             
%             hold on;
%             plot(datenum(2010:0.1:2014.9,0,0), wmm_2015_ngdc(:,3),'r.-','MarkerSize',15);
%             plot(datenum(2010:0.1:2014.9,0,0), wmm_2015_bgs(:,3),'b.-','MarkerSize',15);
%             
%             
%             set(gca,'FontSize',16);
%             ylabel('Z nT');
%             axis([datenum(2010,1,1), datenum(2015,1,1), -inf inf]);
%             datetick('x','keeplimits');
%             
%             
%             set(f,'PaperPositionMode','auto');
%             saveas(f,['/Users/manojnair/projects/wmm2015_validation/' sprintf('%s',spline_obs(i,:)) '_Sec_Var_NGDC_vs_NGDC11'],'epsc');
            
            
            %pause;
            
            clear  sp_x sp_y sp_z;
          %  clf;
            
        end;
        
    end;
    
end;


%% plot maps

% subplot(311);
% worldmap('world');
% plotm(lat,long)
% scatterm(delta_sec(:,1),delta_sec(:,2),40,delta_sec(:,3),'filled');
% caxis([-15,15]);
% title('X OBS-WMM2015 (nT/Yr)')
% 
% colorbar;
% subplot(312)
% worldmap('world');
% plotm(lat,long)
% scatterm(delta_sec(:,1),delta_sec(:,2),40,delta_sec(:,4),'filled');
% caxis([-15,15]);
% title('Y OBS-WMM2015 (nT/Yr)')
% 
% colorbar
% subplot(313)
% worldmap('world');
% plotm(lat,long)
% scatterm(delta_sec(:,1),delta_sec(:,2),40,delta_sec(:,5),'filled');
% caxis([-15,15]);
% title('Z OBS-WMM2015 (nT/Yr)')
% 
% colorbar;
% 
